Processing Pipeline for Pancreatic Endocrinogenesis Data from Bastidas-Ponce et al (2019)

Bacher Group

Author
Affiliation

Jack R. Leary

University of Florida

Published

October 7, 2024

1 Libraries

Before we do anything else we load in the necessary software packages in both R & Python.

1.1 R

Code 
library(dplyr)       # data manipulation
library(GENIE3)      # trajectory GRN construction
library(Seurat)      # scRNA-seq tools
library(scLANE)      # trajectory DE
library(ggplot2)     # pretty plots
library(biomaRt)     # gene annotation
library(tradeSeq)    # more trajectory DE
library(patchwork)   # plot alignment
library(reticulate)  # Python interface
select <- dplyr::select
rename <- dplyr::rename

1.2 Python

Code 
import warnings                                                                 # filter out warnings
import numpy as np                                                              # linear algebra tools
import scanpy as sc                                                             # scRNA-seq tools
import pandas as pd                                                             # DataFrames
import anndata as ad                                                            # scRNA-seq data structures
import scvelo as scv                                                            # RNA velocity
import cellrank as cr                                                           # cell fate estimation
import matplotlib.pyplot as plt                                                 # pretty plots
from scipy.io import mmread, mmwrite                                            # sparse matrix IO
from cellrank.estimators import GPCCA                                           # CellRank estimator
from scipy.sparse import coo_matrix, csr_matrix                                 # sparse matrices
from cellrank.kernels import PseudotimeKernel, CytoTRACEKernel, VelocityKernel  # CellRank kernels

Here we set a few global variables to reduce the amount of printed output caused by our Python code.

Code 
warnings.simplefilter('ignore', category=UserWarning)
sc.settings.verbosity = 0
scv.settings.verbosity = 0
cr.settings.verbosity = 0

2 Visualization tools

To make our visualizations prettier we’ll define some consistent color palettes, and set a theme for matplotlib that matches ggplot2::theme_classic().

2.1 Color palettes

Code 
palette_heatmap <- paletteer::paletteer_d("wesanderson::Zissou1")
palette_cluster <- paletteer::paletteer_d("ggsci::category20_d3")
palette_celltype <- c("#A82203FF", "#208CC0FF", "#F1AF3AFF", "#CF5E4EFF", "#00991AFF", "#003967FF", "#6BD76BFF", "#660099FF")
names(palette_celltype) <- c("Ductal", "Ngn3 low EP", "Ngn3 high EP", "Pre-endocrine", "Beta", "Alpha", "Delta", "Epsilon")

2.2 Theme for matplotlib

Code 
base_size = 12
plt.rcParams.update({
    # font
    'font.size': base_size, 
    'font.weight': 'normal',
    # figure
    'figure.dpi': 320, 
    'figure.edgecolor': 'white', 
    'figure.facecolor': 'white', 
    'figure.figsize': (6, 4), 
    'figure.constrained_layout.use': True,
    # axes
    'axes.edgecolor': 'black',
    'axes.grid': False,
    'axes.labelpad': 2.75,
    'axes.labelsize': base_size * 0.8,
    'axes.linewidth': 1.5,
    'axes.spines.right': False,
    'axes.spines.top': False,
    'axes.titlelocation': 'left',
    'axes.titlepad': 11,
    'axes.titlesize': base_size,
    'axes.titleweight': 'normal',
    'axes.xmargin': 0.1, 
    'axes.ymargin': 0.1, 
    # legend
    'legend.borderaxespad': 1,
    'legend.borderpad': 0.5,
    'legend.columnspacing': 2,
    'legend.fontsize': base_size * 0.8,
    'legend.frameon': False,
    'legend.handleheight': 1,
    'legend.handlelength': 1.2,
    'legend.labelspacing': 1,
    'legend.title_fontsize': base_size, 
    'legend.markerscale': 1.5
})

3 Helper functions

We first define a utility function to make our plot legends cleaner to read & easier to make.

Code 
guide_umap <- function(key.size = 4) {
  ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(size = key.size,
                                                                    alpha = 1, 
                                                                    stroke = 0.25)))
}

4 Data

4.1 Table of murine transcription factors

First we read in a table of Mus musculus transcription factors (TFs) from Lambert et al (2018). Using the biomaRt package we connect to Ensembl and add

Code 
mm_ensembl <- useMart("ensembl", 
                      dataset = "mmusculus_gene_ensembl", 
                      host = "https://useast.ensembl.org")
mm_tf_raw <- readr::read_tsv("https://raw.githubusercontent.com/gifford-lab/ReprogrammingRecovery/main/data/mouse_ensemble_tfs_from_lambertetal_isyes.unique.txt",
                             num_threads = 1,
                             col_names = FALSE,
                             show_col_types = FALSE) %>%
             magrittr::set_colnames("ensembl_id")
mm_tfs <- getBM(attributes = c("ensembl_gene_id", "mgi_symbol", "entrezgene_id", "description", "gene_biotype"),
                filters = "ensembl_gene_id",
                values = mm_tf_raw$ensembl_id,
                mart = mm_ensembl,
                uniqueRows = TRUE) %>%
          rename(ensembl_id = ensembl_gene_id,
                 entrez_id = entrezgene_id) %>%
          arrange(ensembl_id) %>%
          mutate(mgi_symbol = if_else(mgi_symbol == "", NA_character_, mgi_symbol),
                 description = gsub("\\[Source:.*", "", description))

4.2 Pancreatic endocrinogenesis dataset

We load the well-known day E15.5 pancreatic endocrinogenesis data from Bastidas-Ponce et al (2019), which is often used as a good benchmark dataset due to the simplicity of the underlying trajectory manifold.

Code 
ad_panc = scv.datasets.pancreas()
Code 
ad_panc.obs.rename(columns={'clusters': 'celltype', 'clusters_coarse': 'celltype_coarse'}, inplace=True)
ad_panc.uns['celltype_colors'] = np.array(['#A82203FF', '#208CC0FF', '#F1AF3AFF', '#CF5E4EFF', '#00991AFF', '#003967FF', '#6BD76BFF', '#660099FF'])
ad_panc.layers['spliced_counts'] = ad_panc.layers['spliced'].copy()
ad_panc.layers['unspliced_counts'] = ad_panc.layers['unspliced'].copy()
ad_panc.raw = ad_panc

5 Preprocessing

We’ll begin our analysis by running our cells through a standard preprocessing pipeline composed of QC, normalization, HVG selection, dimension reduction, & clustering.

5.1 QC

We filter out cells with a (spliced) sequencing depth of less than 1,000, then remove out genes expressed in less than 5 cells.

Code 
sc.pp.filter_cells(ad_panc, min_counts=1000)
sc.pp.filter_genes(ad_panc, min_cells=5)

5.2 HVG selection

Here we select the top 3000 most highly-variable genes (HVGs), which we’ll use as input to PCA.

Code 
sc.pp.highly_variable_genes(
    ad_panc, 
    n_top_genes=3000, 
    flavor='seurat_v3', 
    subset=False
)

5.3 Normalization

We next depth- and log1p-normalize the spliced mRNA counts matrix.

Code 
ad_panc.X = sc.pp.normalize_total(ad_panc, target_sum=1e4, inplace=False)['X']
sc.pp.log1p(ad_panc)

5.4 PCA embedding

Using PCA we generate a 50-dimensional linear reduction of the normalized spliced counts.

Code 
sc.pp.scale(ad_panc)
sc.tl.pca(
    ad_panc, 
    n_comps=50, 
    random_state=312, 
    use_highly_variable=True
)

Plotting the first two PCs shows clear separation by (coarse) celltype with e.g., variation in the mature endocrine celltypes being less clear.

Code 
sc.pl.embedding(
    ad_panc, 
    basis='pca', 
    color='celltype',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.75,
    show=False
)
plt.gca().set_xlabel('PC 1')
plt.gca().set_ylabel('PC 2')
plt.show()
Figure 1: PCA embedding colored by celltype.

5.5 SNN graph estimation

Next, we identify the 20 nearest-neighbors (NNs) for each cell using the cosine distance in PCA space.

Code 
sc.pp.neighbors(
    ad_panc, 
    n_neighbors=20,
    n_pcs=30,  
    metric='cosine', 
    random_state=312
)

5.6 Clustering

Using the Leiden algorithm we partition the SNN graph into clusters.

Code 
sc.tl.leiden(
    ad_panc, 
    resolution=0.3, 
    random_state=312
)

The identified clusters roughly match our celltypes.

Code 
sc.pl.embedding(
    ad_panc, 
    basis='pca', 
    color='leiden',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.75,
    show=False
)
plt.gca().set_xlabel('PC 1')
plt.gca().set_ylabel('PC 2')
plt.show()
Figure 2: PCA embedding colored by Leiden cluster.

5.7 Moment-based imputation

Moving on, we create smoothed versions of the spliced & unspliced counts across each cell’s 20 NNs.

Code 
scv.pp.moments(
    ad_panc, 
    n_pcs=None,
    n_neighbors=20
)

5.8 Undirected PAGA embedding

We use the PAGA algorithm to generate a graph abstraction of the relationships between celltypes.

Code 
sc.tl.paga(ad_panc, groups='celltype')

Visualizing the results shows a branching structure that matches what we expect biologically.

Code 
sc.pl.paga(
    ad_panc, 
    fontoutline=True, 
    fontsize=12, 
    frameon=True, 
    random_state=312, 
    show=False
)
plt.gca().set_xlabel('PAGA 1')
plt.gca().set_ylabel('PAGA 2')
plt.show()
Figure 3: PAGA embedding colored by celltype. Lines are weighted by (undirected) connectivity strength.

5.9 UMAP embedding

We then create a nonlinear dimension reduction of the data using the UMAP algorithm.

Code 
sc.tl.umap(
    ad_panc, 
    init_pos='paga', 
    random_state=312
)

In UMAP space the transitions between celltypes are clear and well-represented.

Code 
sc.pl.embedding(
    ad_panc, 
    basis='umap', 
    color='celltype',
    title='',
    frameon=True, 
    size=30, 
    alpha=0.75,
    show=False
)
plt.gca().set_xlabel('UMAP 1')
plt.gca().set_ylabel('UMAP 2')
plt.show()
Figure 4: UMAP embedding colored by celltype.

5.10 Force-directed graph embedding

Using the ForceAtlas2 algorithm we generate a 2-dimensional graph layout of our cells.

Code 
sc.tl.draw_graph(
    ad_panc, 
    layout='fa',
    init_pos='paga', 
    random_state=312
)

While the force-directed graph embedding does recapitulate the differences and transitions between celltypes, it’s a little noisier than the UMAP embedding, so we’ll stick with UMAP instead going forward.

Code 
sc.pl.draw_graph(
    ad_panc, 
    color='celltype', 
    title='', 
    alpha=0.75, 
    size=30, 
    show=False, 
    frameon=True
)
plt.gca().set_xlabel('FA 1')
plt.gca().set_ylabel('FA 2')
plt.show()
Figure 5: Force-directed graph embedding colored by celltype.

5.11 Diffusion map embedding

Our last embedding will be generated using diffusion maps, which explicitly attempts to preserve transitional structures in the data.

Code 
sc.tl.diffmap(
    ad_panc, 
    random_state=312, 
    n_comps=16
)
ad_panc.obsm['X_diffmap_old'] = ad_panc.obsm['X_diffmap']
ad_panc.obsm['X_diffmap'] = ad_panc.obsm['X_diffmap'][:, 1:]

While the diffusion map embedding captures some characteristics of the data, it masks variation in the mature endocrine celltypes.

Code 
sc.pl.embedding(
    ad_panc, 
    basis='diffmap', 
    color='celltype',
    title='', 
    frameon=True, 
    alpha=0.75, 
    size=30, 
    components='1,2', 
    show=False
)
plt.gca().set_xlabel('DC 1')
plt.gca().set_ylabel('DC 2')
plt.show()
Figure 6: Diffusion map embedding colored by celltype.

5.12 Diffusion pseudotime estimation

We next compute a diffusion pseudotime estimate for each cell using the first diffusion component.

Code 
ad_panc.uns['iroot'] = np.argmin(ad_panc.obsm['X_diffmap'][:, 0])
sc.tl.dpt(ad_panc, n_dcs=1)

Overall the diffusion pseudotime captures endocrinogenesis progression well, but it assigns very similar values to all the mature endocrine celltypes which isn’t ideal.

Code 
sc.pl.embedding(
    ad_panc, 
    basis='diffmap', 
    color='dpt_pseudotime',
    title='', 
    frameon=True, 
    alpha=0.75, 
    size=30, 
    components='1,2', 
    show=False
)
plt.gca().set_xlabel('DC 1')
plt.gca().set_ylabel('DC 2')
plt.show()
Figure 7: Diffusion map embedding colored by diffusion pseudotime.

6 Analysis

We’ll perform several modes of analysis, namely RNA velocity analysis with scVelo, trajectory DE testing with scLANE, and cell fate analysis with CellRank.

6.1 RNA velocity estimation

We run the dynamic velocity model, estimate a velocity graph, and embed the velocity vectors in UMAP space. We’ll also record the runtime in order to compare it with that of scLANE.

Code 
velo_start <- Sys.time()
Code 
scv.tl.recover_dynamics(ad_panc, n_jobs=8)
scv.tl.velocity(
    ad_panc, 
    mode='dynamical', 
    use_highly_variable=True, 
    filter_genes=False
)
scv.tl.velocity_graph(
  ad_panc, 
  compute_uncertainties=True, 
  n_neighbors=20, 
  n_jobs=8
)
scv.tl.velocity_embedding(ad_panc, basis='umap')

We then estimate the confidence of the velocities along with a gene-shared latent time ordering.

Code 
scv.tl.velocity_confidence(ad_panc)
scv.tl.latent_time(ad_panc)

The total runtime was:

Code 
velo_end <- Sys.time()
velo_diff <- velo_end - velo_start
velo_diff
Time difference of 2.758425 mins

The latent time representation seems to capture well the progression from endocrine progenitor cells to mature endocrine phenotypes, with beta cells having the highest latent time measurements overall.

Code 
sc.pl.embedding(
    ad_panc, 
    basis='umap', 
    color='latent_time',
    title='', 
    frameon=True, 
    alpha=0.75, 
    size=30,  
    show=False
)
plt.gca().set_xlabel('UMAP 1')
plt.gca().set_ylabel('UMAP 2')
plt.show()
Figure 8: UMAP embedding colored by latent time.

The velocity streamlines show clear progression through the endocrine progenitors up through the mature alpha, beta, delta, and epsilon celltypes. However, directionality is less clear in the ductal cells, most likely due to a combination of cell cycle effects and our explicit incorporation of basal transcription in the velocity model.

Code 
scv.pl.velocity_embedding_stream(
    ad_panc, 
    basis='umap', 
    color='celltype',
    alpha=0.75, 
    size=30, 
    frameon=True, 
    linewidth=1, 
    legend_loc='right margin',
    arrow_size=1,
    title='', 
    show=False
)
plt.gca().set_xlabel('UMAP 1')
plt.gca().set_ylabel('UMAP 2')
plt.show()
Figure 9: UMAP embedding colored by celltype overlaid with velocity vector streamlines.

6.1.1 Directed PAGA embedding

Here we use the PAGA algorithm to generate a directed graph representation of the relationships between celltypes using RNA velocity as prior information.

Code 
scv.tl.paga(
  ad_panc, 
  groups='celltype', 
  use_time_prior='latent_time'
)

The RNA velocity measurements have allowed us to identify directed relationships between celltypes, which we visualize below.

Code 
scv.pl.paga(
    ad_panc, 
    basis='umap', 
    color='celltype',
    title='',
    frameon=True, 
    min_edge_width=1, 
    random_state=312,
    show=False
)
plt.gca().set_xlabel('UMAP 1')
plt.gca().set_ylabel('UMAP 2')
plt.show()
Figure 10: Directed PAGA embedding of the relationships between celltypes. Linewidth denotes strength of connectivity.

6.2 Save preprocessed data

We save the preprocessed object to an .h5ad file.

Code 
ad_panc.write_h5ad('../../data/pancreas_E15.5/ad_panc.h5ad')

We also save each of the individual pieces of the AnnData object to files, which will allow us to read the data into R & create a Seurat object.

Code 
ad_panc.obs.reset_index(inplace=False).rename(columns={'index': 'cell'}).to_csv('../../data/pancreas_E15.5/conversion_files/obs.csv')
ad_panc.var.reset_index(inplace=False).rename(columns={'index': 'gene'}).to_csv('../../data/pancreas_E15.5/conversion_files/var.csv')
pd.DataFrame(ad_panc.obsm['X_diffmap']).reset_index(inplace=False).to_csv('../../data/pancreas_E15.5/conversion_files/diffmap.csv')
pd.DataFrame(ad_panc.obsm['X_pca']).reset_index(inplace=False).to_csv('../../data/pancreas_E15.5/conversion_files/PCA.csv')
pd.DataFrame(ad_panc.obsm['X_draw_graph_fa']).reset_index(inplace=False).to_csv('../../data/pancreas_E15.5/conversion_files/FA.csv')
pd.DataFrame(ad_panc.obsm['X_umap']).reset_index(inplace=False).to_csv('../../data/pancreas_E15.5/conversion_files/UMAP.csv')
mmwrite('../../data/pancreas_E15.5/conversion_files/spliced_counts.mtx', a=ad_panc.layers['spliced_counts'])
mmwrite('../../data/pancreas_E15.5/conversion_files/unspliced_counts.mtx', a=ad_panc.layers['unspliced_counts'])

6.3 Import data from Python

First we read in the data that we processed with scanpy. Expand the code block for details.

Code 
cell_metadata <- readr::read_csv("../../data/pancreas_E15.5/conversion_files/obs.csv", 
                                 show_col_types = FALSE, 
                                 col_select = -1) %>% 
                 as.data.frame() %>% 
                 magrittr::set_rownames(.$cell)
gene_metadata <- readr::read_csv("../../data/pancreas_E15.5/conversion_files/var.csv", 
                                 show_col_types = FALSE, 
                                 col_select = -1)
embed_PCA <- readr::read_csv("../../data/pancreas_E15.5/conversion_files/PCA.csv",
                             show_col_types = FALSE, 
                             col_select = -c(1:2)) %>% 
             as.data.frame() %>% 
             magrittr::set_colnames(paste0("PC", 1:ncol(.))) %>% 
             magrittr::set_rownames(cell_metadata$cell)
embed_PCA <- CreateDimReducObject(embeddings = as.matrix(embed_PCA),
                                  key = "PC_", 
                                  global = TRUE, 
                                  assay = "RNA")
embed_FA <- readr::read_csv("../../data/pancreas_E15.5/conversion_files/FA.csv", 
                            show_col_types = FALSE, 
                            col_select = -c(1:2)) %>% 
            as.data.frame() %>% 
            magrittr::set_colnames(paste0("FA", 1:ncol(.))) %>% 
            magrittr::set_rownames(cell_metadata$cell)
embed_FA <- CreateDimReducObject(embeddings = as.matrix(embed_FA),
                                 key = "FA_", 
                                 global = TRUE, 
                                 assay = "RNA")
embed_UMAP <- readr::read_csv("../../data/pancreas_E15.5/conversion_files/UMAP.csv", 
                              show_col_types = FALSE, 
                               col_select = -c(1:2)) %>% 
              as.data.frame() %>% 
              magrittr::set_colnames(paste0("UMAP", 1:ncol(.))) %>% 
              magrittr::set_rownames(cell_metadata$cell)
embed_UMAP <- CreateDimReducObject(embeddings = as.matrix(embed_UMAP),
                                   key = "UMAP_", 
                                   global = TRUE, 
                                   assay = "RNA")
embed_diffmap <- readr::read_csv("../../data/pancreas_E15.5/conversion_files/diffmap.csv", 
                                 show_col_types = FALSE, 
                                  col_select = -c(1:2)) %>% 
                 as.data.frame() %>% 
                 magrittr::set_colnames(paste0("DC", 1:ncol(.))) %>% 
                 magrittr::set_rownames(cell_metadata$cell)
embed_diffmap <- CreateDimReducObject(embeddings = as.matrix(embed_diffmap),
                                      key = "DC_", 
                                      global = TRUE, 
                                      assay = "RNA")
spliced_counts <- Matrix::t(Matrix::readMM("../../data/pancreas_E15.5/conversion_files/spliced_counts.mtx"))
unspliced_counts <- Matrix::t(Matrix::readMM("../../data/pancreas_E15.5/conversion_files/unspliced_counts.mtx"))
colnames(spliced_counts) <- colnames(unspliced_counts) <-cell_metadata$cell
rownames(spliced_counts) <- rownames(unspliced_counts) <- gene_metadata$gene
spliced_counts <- CreateAssayObject(counts = spliced_counts,
                                    min.cells = 0,
                                    min.features = 0, 
                                    key = "spliced")
unspliced_counts <- CreateAssayObject(counts = unspliced_counts, 
                                      min.cells = 0, 
                                      min.features = 0, 
                                      key = "unspliced")

Now we can create a Seurat object. After adding both assays and the embeddings, we normalize both the spliced & unspliced counts, then re-add the HVG information (mean, variance, etc.) for each gene.

Code 
seu_panc <- CreateSeuratObject(spliced_counts, 
                               assay = "spliced",
                               meta.data = cell_metadata, 
                               project = "pancreas_E15.5", 
                               min.cells = 0, 
                               min.features = 0)
seu_panc@meta.data <- mutate(seu_panc@meta.data, 
                             celltype = factor(celltype, levels = c("Ductal",
                                                                    "Ngn3 low EP", 
                                                                    "Ngn3 high EP", 
                                                                    "Pre-endocrine", 
                                                                    "Beta", 
                                                                    "Alpha", 
                                                                    "Delta", 
                                                                    "Epsilon")))
seu_panc@assays$unspliced <- unspliced_counts
seu_panc@reductions$pca <- embed_PCA
seu_panc@reductions$diffmap <- embed_diffmap
seu_panc@reductions$fa <- embed_FA
seu_panc@reductions$umap <- embed_UMAP
seu_panc <- seu_panc %>% 
            NormalizeData(assay = "spliced", verbose = FALSE) %>% 
            NormalizeData(assays = "unspliced", verbose = FALSE) %>% 
            FindVariableFeatures(assay = "spliced", 
                                 selection.method = "vst", 
                                 verbose = FALSE)

We can see that our UMAP embedding has been accurately preserved, even down to the colors of the celltypes.

Code 
p0 <- as.data.frame(Embeddings(seu_panc, "umap")) %>% 
      magrittr::set_colnames(c("UMAP_1", "UMAP_2")) %>% 
      mutate(celltype = seu_panc$celltype) %>% 
      ggplot(aes(x = UMAP_1, y = UMAP_2, color = celltype)) + 
      geom_point(size = 1.25, 
                 alpha = 0.75, 
                 stroke = 0) + 
      scale_color_manual(values = palette_celltype) + 
      labs(x = "UMAP 1",  y = "UMAP 2") + 
      theme_scLANE(umap = TRUE) + 
      theme(plot.title = element_blank(), 
            legend.title = element_blank()) + 
      guide_umap()
p0
Figure 11: Preserved UMAP embedding from Python colored by celltype.

6.4 Trajectory DE testing with scLANE

We use the scLANE package to perform trajectory differential expression testing across latent time for the top 3000 candidate genes.

Code 
cell_offset <- createCellOffset(seu_panc)
pt_df <- data.frame(LT = seu_panc$latent_time)
candidate_genes <- chooseCandidateGenes(seu_panc, 
                                        group.by.subject = FALSE,
                                        n.desired.genes = 4000L)
scLANE_models <- testDynamic(seu_panc, 
                             pt = pt_df, 
                             genes = candidate_genes, 
                             size.factor.offset = cell_offset, 
                             n.cores = 16L, 
                             verbose = FALSE)
scLANE testing in GLM mode completed for 4000 genes across 1 lineage in 34.907 mins

We generate a tidy table of DE test statistics, then identify a set of dynamic genes.

Code 
scLANE_de_res_tidy <- getResultsDE(scLANE_models)
dyn_genes <- filter(scLANE_de_res_tidy, Gene_Dynamic_Overall == 1) %>% 
             distinct(Gene) %>% 
             pull(Gene)

6.5 tradeSeq DE analysis

Code 
ts_start <- Sys.time()
bioc_par <- BiocParallel::MulticoreParam(workers = 16L, RNGseed = 312)
spliced_counts <- as.matrix(seu_panc@assays$spliced$counts)[candidate_genes, ]
k_eval <- evaluateK(spliced_counts, 
                    pseudotime = pt_df, 
                    cellWeights = matrix(rep(1, nrow(pt_df)), ncol = 1), 
                    offset = log(1 / cell_offset), 
                    k = 3:10, 
                    plot = FALSE, 
                    nGenes = 500, 
                    verbose = FALSE, 
                    parallel = TRUE, 
                    BPPARAM = bioc_par)
best_k <- c(3:10)[which.min(abs(colMeans(k_eval - rowMeans(k_eval))))]  # choose k w/ lowest MAD from mean AIC 
ts_models <- fitGAM(spliced_counts, 
                    pseudotime = pt_df, 
                    cellWeights = matrix(rep(1, nrow(pt_df)), ncol = 1), 
                    offset = log(1 / cell_offset), 
                    nknots = best_k, 
                    sce = FALSE, 
                    parallel = TRUE, 
                    verbose = FALSE, 
                    BPPARAM = bioc_par)
BiocParallel::bpstop(bioc_par)
ts_de_res <- associationTest(ts_models, global = TRUE) %>% 
                             arrange(desc(waldStat)) %>% 
                             mutate(gene = rownames(.), 
                                    pvalue_adj = p.adjust(pvalue, method = "fdr"), 
                                    gene_dynamic_overall = if_else(pvalue_adj < 0.01, 1, 0)) %>% 
                             relocate(gene)
ts_end <- Sys.time()
ts_diff <- ts_end - ts_start

The runtime for tradeSeq with 5 knots per gene is:

Code 
ts_diff
Time difference of 13.99666 mins

6.6 Cell fate analysis with CellRank

We’ll now use the CellRank package to analyze initial and terminal states in our data.

6.6.1 CytoTRACE kernel

We begin by computing a CytoTRACE kernel; CytoTRACE uses the number of genes expressed in a cell as a proxy measurement for differentiation potential, with the assumption being that cells expressing more genes are more likely to be progenitor celltypes & vice versa.

Code 
ctk = CytoTRACEKernel(ad_panc).compute_cytotrace()

We plot the CytoTRACE scores below, and see that the ductal & endocrine progenitor populations score most highly (as expected).

Code 
sc.pl.embedding(
    ad_panc,
    color=['ct_score', 'ct_pseudotime'],
    basis='umap',
    color_map='magma', 
    show=False, 
    size=30, 
    alpha=0.75, 
    title=['CytoTRACE score', 'CytoTRACE Pseudotime']
)
plt.show()
Figure 12: UMAP embedding colored by CytoTRACE score and pseudotime. Higher scores indicate higher differentiation potential.

Next, we estimate a cell-cell transition probability matrix.

Code 
ctk.compute_transition_matrix(threshold_scheme='soft', show_progress_bar=False)

Projecting that matrix onto our UMAP embedding reveals the differentiation directions of our cells based on their differentiation potential. The directionality in the mature celltypes is uncertain and in some cases outright wrong, which makes sense as their CytoTRACE scores are all very similar. This can be ameliorated by combining the CytoTRACE score kernel with kernels derived from other modalities that better capture the heterogeneity of the mature endocrine cells.

Code 
ctk.plot_projection(
  basis='umap', 
  recompute=True, 
  color='celltype',
  title='',
  legend_loc='right margin', 
  size=30, 
  alpha=0.75, 
  linewidth=2, 
  frameon=True, 
  show=False
)
plt.gca().set_xlabel('UMAP 1')
plt.gca().set_ylabel('UMAP 2')
plt.show()
Figure 13: Streamline embedding plot showing differentiation directions as inferred from CytoTRACE scores.

6.6.2 Velocity kernel

Moving on, we use the RNA velocity estimates to create a velocity kernel. Using Monte Carlo sampling we can propagate forward the velocity uncertainties we computed earlier to the cell-cell transition probability matrix.

Code 
vk = cr.kernels.VelocityKernel(ad_panc)
vk.compute_transition_matrix(
    model='monte_carlo', 
    similarity='cosine', 
    seed=312, 
    show_progress_bar=False
)

The velocity kernel projection of course looks highly similar to the velocity streamline embedding computed after running scVelo.

Code 
vk.plot_projection(
  basis='umap', 
  recompute=True, 
  color='celltype',
  title='',
  legend_loc='right margin', 
  size=30, 
  alpha=0.75, 
  linewidth=2, 
  frameon=True, 
  show=False
)
plt.gca().set_xlabel('UMAP 1')
plt.gca().set_ylabel('UMAP 2')
plt.show()
Figure 14: Streamline embedding plot showing differentiation directions as inferred from RNA velocity vectors.

6.6.3 Pseudotime kernel

Lastly, we create a pseudotime-based kernel using diffusion pseudotime, then compute another cell-cell transition probability matrix.

Code 
pk = PseudotimeKernel(ad_panc, time_key='dpt_pseudotime')
pk.compute_transition_matrix(threshold_scheme='soft', show_progress_bar=False)

Plotting the matrix’s projection onto our UMAP embedding shows smooth transitions from the ductal cells through the endocrine progenitors towards the mature endocrine celltypes.

Code 
pk.plot_projection(
  basis='umap', 
  recompute=True, 
  color='celltype',
  title='',
  legend_loc='right margin', 
  size=30, 
  alpha=0.75, 
  linewidth=2, 
  frameon=True, 
  show=False
)
plt.gca().set_xlabel('UMAP 1')
plt.gca().set_ylabel('UMAP 2')
plt.show()
Figure 15: UMAP embedding overlaid with streamline arrows derived from diffusion pseudotime.

6.6.4 Combined kernel

A key feature of CellRank is the ability to combine kernels in a weighted manner, which we do so below using unequal weights. This provides us with a coherent cell-cell transition probability matrix derived from multiple data modalities.

Code 
ck = 0.25 * ctk + 0.5 * vk + 0.25 * pk

We visualize the combined projection below. It’s evident that combining the different modalities has resulted in a harmonized, accurate representation of the transitions between celltypes.

Code 
ck.plot_projection(
  basis='umap', 
  recompute=True, 
  color='celltype',
  title='',
  legend_loc='right margin', 
  size=30, 
  alpha=0.75, 
  linewidth=2, 
  frameon=True, 
  show=False
)
plt.gca().set_xlabel('UMAP 1')
plt.gca().set_ylabel('UMAP 2')
plt.show()
Figure 16: UMAP embedding overlaid with streamline arrows derived from the combined kernel.

We can also simulate random walks on the high-dimensional cell-cell graph starting from the ductal population. We see that the random walks tend to terminate in the mature endocrine celltypes, as expected.

Code 
ck.plot_random_walks(
    n_sims=100,
    start_ixs={'celltype': 'Ductal'},
    basis='umap',
    color='celltype',
    legend_loc='right margin',
    seed=312, 
    frameon=True, 
    title='', 
    linewidth=0.5, 
    linealpha=0.25, 
    ixs_legend_loc='upper', 
    show_progress_bar=False
)
plt.gca().set_xlabel('UMAP 1')
plt.gca().set_ylabel('UMAP 2')
plt.show()
Figure 17: UMAP embedding overlaid with random walks along the cell-cell graph derived from the combined kernel. Starting points are highlighted in black and ending points in yellow.

7 Fitting a gene regulatory network

Finally, in order to validate the links between transcription factors and target genes, we build a gene regulatory network (GRN) using the GENIE3 package. This process takes a while to complete, so we utilize 24 threads to speed things up.

Code 
spliced_norm_counts <- as.matrix(seu_panc@assays$spliced$data[candidate_genes, ])
valid_tfs <- mm_tfs$mgi_symbol[mm_tfs$mgi_symbol %in% rownames(spliced_norm_counts)]
valid_targets <- rownames(spliced_norm_counts)[!rownames(spliced_norm_counts) %in% mm_tfs$mgi_symbol]
grn <- GENIE3(spliced_norm_counts, 
              regulators = valid_tfs, 
              targets = valid_targets, 
              treeMethod = "RF", 
              K = "sqrt", 
              nCores = 36L,
              verbose = FALSE)
grn_df <- getLinkList(grn)

A sample of the GRN results shows us links between TFs and target genes:

Code 
slice_sample(grn_df, n = 10) %>% 
  kableExtra::kable(digits = 4, 
                    row.names = FALSE, 
                    col.names = c("TF", "Target", "Weight"), 
                    booktabs = TRUE) %>% 
  kableExtra::kable_classic(full_width = FALSE, "hover")
TF Target Weight
Foxo1 Ctnnbl1 0.0036
Tet3 Sdhaf4 0.0042
Neurod2 Fam134a 0.0030
Zfp24 Prom1 0.0032
Zfp414 Naa50 0.0045
Mbd3 Hsd17b12 0.0090
Nr3c1 Trim8 0.0073
Irf3 Cmtm8 0.0022
Mbd3 Psmg1 0.0079
Rest Pla2g2f 0.0095
Table 1: A random sample of 10 links between TFs and target genes as estimated by GENIE3.

8 Save data

Code 
saveRDS(seu_panc, file = "../../data/pancreas_E15.5/seu_panc.Rds")
saveRDS(scLANE_models, file = "../../data/pancreas_E15.5/scLANE_models.Rds")
saveRDS(ts_models, file = "../../data/pancreas_E15.5/ts_models.Rds")
saveRDS(grn, file = "../../data/pancreas_E15.5/grn.Rds")

9 Session info

Code 
sessioninfo::session_info()
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.3 (2024-02-29)
 os       Ubuntu 22.04.5 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/New_York
 date     2024-10-07
 pandoc   2.9.2.1 @ /usr/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package              * version    date (UTC) lib source
 abind                  1.4-8      2024-09-12 [2] CRAN (R 4.3.3)
 AnnotationDbi          1.62.2     2023-07-02 [2] Bioconductor
 backports              1.5.0      2024-05-23 [1] CRAN (R 4.2.3)
 bigassertr             0.1.6      2023-01-10 [1] CRAN (R 4.3.2)
 bigparallelr           0.3.2      2021-10-02 [1] CRAN (R 4.3.2)
 bigstatsr              1.5.12     2022-10-14 [1] CRAN (R 4.3.2)
 Biobase                2.60.0     2023-04-25 [2] Bioconductor
 BiocFileCache          2.13.1     2024-10-05 [1] Github (Bioconductor/BiocFileCache@30abbe4)
 BiocGenerics           0.46.0     2023-04-25 [2] Bioconductor
 BiocParallel           1.34.2     2023-05-22 [2] Bioconductor
 biomaRt              * 2.56.1     2023-06-09 [1] Bioconductor
 Biostrings             2.68.1     2023-05-16 [2] Bioconductor
 bit                    4.5.0      2024-09-20 [2] CRAN (R 4.3.3)
 bit64                  4.5.2      2024-09-22 [2] CRAN (R 4.3.3)
 bitops                 1.0-8      2024-07-29 [2] CRAN (R 4.3.3)
 blob                   1.2.4      2023-03-17 [2] CRAN (R 4.3.1)
 boot                   1.3-31     2024-08-28 [4] CRAN (R 4.3.3)
 broom                  1.0.6      2024-05-17 [1] CRAN (R 4.2.3)
 broom.mixed            0.2.9.5    2024-04-01 [1] CRAN (R 4.2.3)
 cachem                 1.1.0      2024-05-16 [2] CRAN (R 4.3.3)
 cli                    3.6.2      2023-12-11 [1] CRAN (R 4.2.3)
 cluster                2.1.6      2023-12-01 [4] CRAN (R 4.3.2)
 coda                   0.19-4.1   2024-01-31 [2] CRAN (R 4.3.2)
 codetools              0.2-20     2024-03-31 [4] CRAN (R 4.3.3)
 colorspace             2.1-1      2024-07-26 [2] CRAN (R 4.3.3)
 cowplot                1.1.3      2024-01-22 [2] CRAN (R 4.3.3)
 crayon                 1.5.3      2024-06-20 [2] CRAN (R 4.3.3)
 curl                   5.2.3      2024-09-20 [2] CRAN (R 4.3.3)
 data.table             1.16.0     2024-08-27 [2] CRAN (R 4.3.3)
 DBI                    1.2.3      2024-06-02 [2] CRAN (R 4.3.3)
 dbplyr                 2.5.0      2024-03-19 [1] CRAN (R 4.3.3)
 DelayedArray           0.26.7     2023-07-28 [2] Bioconductor
 deldir                 2.0-4      2024-02-28 [2] CRAN (R 4.3.3)
 digest                 0.6.35     2024-03-11 [1] CRAN (R 4.2.3)
 doParallel             1.0.17     2022-02-07 [1] CRAN (R 4.2.3)
 doRNG                * 1.8.6      2023-01-16 [2] CRAN (R 4.3.1)
 doSNOW                 1.0.20     2022-02-04 [1] CRAN (R 4.3.2)
 dotCall64              1.1-1      2023-11-28 [2] CRAN (R 4.3.1)
 dplyr                * 1.1.4      2023-11-17 [1] CRAN (R 4.2.3)
 edgeR                  3.42.4     2023-05-31 [2] Bioconductor
 emmeans                1.10.2     2024-05-20 [1] CRAN (R 4.3.3)
 estimability           1.5.1      2024-05-12 [1] CRAN (R 4.3.3)
 evaluate               0.24.0     2024-06-10 [1] CRAN (R 4.3.3)
 fansi                  1.0.6      2023-12-08 [2] CRAN (R 4.3.1)
 farver                 2.1.2      2024-05-13 [1] CRAN (R 4.2.3)
 fastDummies            1.7.4      2024-08-16 [2] CRAN (R 4.3.3)
 fastmap                1.2.0      2024-05-15 [2] CRAN (R 4.3.3)
 ff                     4.0.12     2024-01-12 [1] CRAN (R 4.3.2)
 filelock               1.0.3      2023-12-11 [2] CRAN (R 4.3.3)
 fitdistrplus           1.2-1      2024-07-12 [2] CRAN (R 4.3.3)
 flock                  0.7        2016-11-12 [1] CRAN (R 4.3.2)
 forcats                1.0.0      2023-01-29 [1] CRAN (R 4.3.2)
 foreach              * 1.5.2      2022-02-02 [1] CRAN (R 4.2.3)
 furrr                  0.3.1      2022-08-15 [1] CRAN (R 4.3.2)
 future                 1.33.2     2024-03-26 [1] CRAN (R 4.2.3)
 future.apply           1.11.2     2024-03-28 [2] CRAN (R 4.3.3)
 gamlss                 5.4-22     2024-03-20 [1] CRAN (R 4.3.2)
 gamlss.data            6.0-6      2024-03-14 [1] CRAN (R 4.3.2)
 gamlss.dist            6.1-1      2023-08-23 [1] CRAN (R 4.3.2)
 geeM                   0.10.1     2018-06-18 [1] CRAN (R 4.3.2)
 generics               0.1.3      2022-07-05 [1] CRAN (R 4.2.0)
 GENIE3               * 1.22.0     2023-04-25 [1] Bioconductor
 GenomeInfoDb           1.36.4     2023-10-02 [2] Bioconductor
 GenomeInfoDbData       1.2.10     2023-07-26 [2] Bioconductor
 GenomicRanges          1.52.1     2023-10-08 [2] Bioconductor
 ggplot2              * 3.5.1      2024-04-23 [1] CRAN (R 4.2.3)
 ggrepel                0.9.2      2022-11-06 [1] CRAN (R 4.2.2)
 ggridges               0.5.6      2024-01-23 [2] CRAN (R 4.3.3)
 glm2                 * 1.2.1      2018-08-11 [1] CRAN (R 4.3.2)
 glmmTMB                1.1.9      2024-03-20 [1] CRAN (R 4.3.3)
 globals                0.16.3     2024-03-08 [2] CRAN (R 4.3.3)
 glue                   1.8.0      2024-09-30 [2] CRAN (R 4.3.3)
 goftest                1.2-3      2021-10-07 [2] CRAN (R 4.3.1)
 gridExtra              2.3        2017-09-09 [2] CRAN (R 4.3.1)
 gtable                 0.3.5      2024-04-22 [2] CRAN (R 4.3.3)
 highr                  0.11       2024-05-26 [2] CRAN (R 4.3.3)
 hms                    1.1.3      2023-03-21 [2] CRAN (R 4.3.1)
 htmltools              0.5.8.1    2024-04-04 [1] CRAN (R 4.2.3)
 htmlwidgets            1.6.4      2023-12-06 [2] CRAN (R 4.3.3)
 httpuv                 1.6.15     2024-03-26 [2] CRAN (R 4.3.3)
 httr                   1.4.7      2023-08-15 [2] CRAN (R 4.3.3)
 ica                    1.0-3      2022-07-08 [2] CRAN (R 4.3.1)
 igraph                 2.0.3      2024-03-13 [1] CRAN (R 4.2.3)
 IRanges                2.34.1     2023-06-22 [2] Bioconductor
 irlba                  2.3.5.1    2022-10-03 [2] CRAN (R 4.3.1)
 iterators              1.0.14     2022-02-05 [2] CRAN (R 4.3.1)
 jsonlite               1.8.8      2023-12-04 [1] CRAN (R 4.2.3)
 kableExtra             1.4.0      2024-01-24 [1] CRAN (R 4.2.3)
 KEGGREST               1.40.1     2023-09-29 [2] Bioconductor
 KernSmooth             2.23-24    2024-05-17 [4] CRAN (R 4.3.3)
 knitr                  1.48       2024-07-07 [2] CRAN (R 4.3.3)
 labeling               0.4.3      2023-08-29 [2] CRAN (R 4.3.1)
 later                  1.3.2      2023-12-06 [2] CRAN (R 4.3.3)
 lattice                0.22-6     2024-03-20 [4] CRAN (R 4.3.3)
 lazyeval               0.2.2      2019-03-15 [2] CRAN (R 4.3.1)
 leiden                 0.4.3.1    2023-11-17 [2] CRAN (R 4.3.1)
 lifecycle              1.0.4      2023-11-07 [2] CRAN (R 4.3.1)
 limma                  3.56.2     2023-06-04 [2] Bioconductor
 listenv                0.9.1      2024-01-29 [2] CRAN (R 4.3.3)
 lme4                   1.1-35.5   2024-07-03 [2] CRAN (R 4.3.3)
 lmtest                 0.9-40     2022-03-21 [2] CRAN (R 4.3.1)
 locfit                 1.5-9.10   2024-06-24 [2] CRAN (R 4.3.3)
 magrittr             * 2.0.3      2022-03-30 [1] CRAN (R 4.2.0)
 MASS                   7.3-60     2023-05-04 [4] CRAN (R 4.3.1)
 Matrix                 1.6-5      2024-01-11 [1] CRAN (R 4.2.3)
 MatrixGenerics         1.12.3     2023-07-30 [2] Bioconductor
 matrixStats            1.4.1      2024-09-08 [2] CRAN (R 4.3.3)
 memoise                2.0.1      2021-11-26 [2] CRAN (R 4.3.1)
 mgcv                   1.9-1      2023-12-21 [4] CRAN (R 4.3.2)
 mime                   0.12       2021-09-28 [2] CRAN (R 4.3.1)
 miniUI                 0.1.1.1    2018-05-18 [2] CRAN (R 4.3.1)
 minqa                  1.2.7      2024-05-20 [1] CRAN (R 4.2.3)
 munsell                0.5.1      2024-04-01 [2] CRAN (R 4.3.3)
 mvtnorm                1.3-1      2024-09-03 [2] CRAN (R 4.3.3)
 nlme                   3.1-166    2024-08-14 [4] CRAN (R 4.3.3)
 nloptr                 2.1.1      2024-06-25 [2] CRAN (R 4.3.3)
 numDeriv               2016.8-1.1 2019-06-06 [2] CRAN (R 4.3.1)
 paletteer              1.6.0      2024-01-21 [1] CRAN (R 4.2.3)
 parallelly             1.37.1     2024-02-29 [1] CRAN (R 4.2.3)
 patchwork            * 1.2.0      2024-01-08 [1] CRAN (R 4.2.3)
 pbapply                1.7-2      2023-06-27 [2] CRAN (R 4.3.1)
 pillar                 1.9.0      2023-03-22 [2] CRAN (R 4.3.1)
 pkgconfig              2.0.3      2019-09-22 [2] CRAN (R 4.3.1)
 plotly                 4.10.4     2024-01-13 [2] CRAN (R 4.3.3)
 plyr                   1.8.9      2023-10-02 [2] CRAN (R 4.3.3)
 png                    0.1-8      2022-11-29 [2] CRAN (R 4.3.1)
 polyclip               1.10-7     2024-07-23 [2] CRAN (R 4.3.3)
 prettyunits            1.2.0      2023-09-24 [2] CRAN (R 4.3.3)
 princurve              2.1.6      2021-01-18 [1] CRAN (R 4.3.2)
 prismatic              1.1.2      2024-04-10 [1] CRAN (R 4.2.3)
 progress               1.2.3      2023-12-06 [2] CRAN (R 4.3.3)
 progressr              0.14.0     2023-08-10 [2] CRAN (R 4.3.2)
 promises               1.3.0      2024-04-05 [2] CRAN (R 4.3.3)
 ps                     1.8.0      2024-09-12 [2] CRAN (R 4.3.3)
 purrr                  1.0.2      2023-08-10 [1] CRAN (R 4.2.3)
 R6                     2.5.1      2021-08-19 [2] CRAN (R 4.3.1)
 RANN                   2.6.2      2024-08-25 [2] CRAN (R 4.3.3)
 rappdirs               0.3.3      2021-01-31 [2] CRAN (R 4.3.1)
 RColorBrewer           1.1-3      2022-04-03 [2] CRAN (R 4.3.1)
 Rcpp                   1.0.13     2024-07-17 [2] CRAN (R 4.3.3)
 RcppAnnoy              0.0.22     2024-01-23 [2] CRAN (R 4.3.3)
 RcppEigen              0.3.4.0.0  2024-02-28 [1] CRAN (R 4.2.3)
 RcppHNSW               0.6.0      2024-02-04 [2] CRAN (R 4.3.3)
 RCurl                  1.98-1.16  2024-07-11 [2] CRAN (R 4.3.3)
 readr                  2.1.5      2024-01-10 [2] CRAN (R 4.3.3)
 rematch2               2.1.2      2020-05-01 [2] CRAN (R 4.3.1)
 reshape2               1.4.4      2020-04-09 [2] CRAN (R 4.3.1)
 reticulate           * 1.35.0     2024-01-31 [1] CRAN (R 4.2.3)
 rlang                  1.1.4      2024-06-04 [1] CRAN (R 4.3.3)
 rmarkdown              2.27       2024-05-17 [1] CRAN (R 4.3.3)
 rmio                   0.4.0      2022-02-17 [1] CRAN (R 4.3.2)
 rngtools             * 1.5.2      2021-09-20 [2] CRAN (R 4.3.1)
 ROCR                   1.0-11     2020-05-02 [2] CRAN (R 4.3.1)
 RSpectra               0.16-2     2024-07-18 [2] CRAN (R 4.3.3)
 RSQLite                2.3.7      2024-05-27 [2] CRAN (R 4.3.3)
 rstudioapi             0.16.0     2024-03-24 [1] CRAN (R 4.2.3)
 Rtsne                  0.17       2023-12-07 [2] CRAN (R 4.3.3)
 S4Arrays               1.2.1      2024-03-04 [1] Bioconductor 3.18 (R 4.3.2)
 S4Vectors              0.38.2     2023-09-22 [2] Bioconductor
 scales                 1.3.0      2023-11-28 [1] CRAN (R 4.2.3)
 scattermore            1.2        2023-06-12 [2] CRAN (R 4.3.1)
 scLANE               * 0.8.4      2024-10-05 [1] Github (jr-leary7/scLANE@d1eb84d)
 sctransform            0.4.1      2023-10-19 [1] CRAN (R 4.2.3)
 sessioninfo            1.2.2      2021-12-06 [2] CRAN (R 4.3.1)
 Seurat               * 5.1.0      2024-05-10 [1] CRAN (R 4.3.3)
 SeuratObject         * 5.0.2      2024-05-08 [1] CRAN (R 4.3.3)
 shiny                  1.9.1      2024-08-01 [2] CRAN (R 4.3.3)
 SingleCellExperiment   1.22.0     2023-04-25 [2] Bioconductor
 slingshot              2.6.0      2022-11-01 [1] Bioconductor
 snow                   0.4-4      2021-10-27 [2] CRAN (R 4.3.1)
 sp                   * 2.1-4      2024-04-30 [3] CRAN (R 4.3.3)
 spam                   2.10-0     2023-10-23 [2] CRAN (R 4.3.1)
 spatstat.data          3.1-2      2024-06-21 [2] CRAN (R 4.3.3)
 spatstat.explore       3.3-2      2024-08-21 [2] CRAN (R 4.3.3)
 spatstat.geom          3.3-3      2024-09-18 [2] CRAN (R 4.3.3)
 spatstat.random        3.3-2      2024-09-18 [2] CRAN (R 4.3.3)
 spatstat.sparse        3.1-0      2024-06-21 [2] CRAN (R 4.3.3)
 spatstat.univar        3.0-1      2024-09-05 [2] CRAN (R 4.3.3)
 spatstat.utils         3.1-0      2024-08-17 [2] CRAN (R 4.3.3)
 stringi                1.8.4      2024-05-06 [2] CRAN (R 4.3.3)
 stringr                1.5.1      2023-11-14 [2] CRAN (R 4.3.2)
 SummarizedExperiment   1.30.2     2023-06-06 [2] Bioconductor
 survival               3.7-0      2024-06-05 [4] CRAN (R 4.3.3)
 svglite                2.1.3      2023-12-08 [1] CRAN (R 4.3.3)
 systemfonts            1.1.0      2024-05-15 [2] CRAN (R 4.3.3)
 tensor                 1.5        2012-05-05 [2] CRAN (R 4.3.1)
 tibble                 3.2.1      2023-03-20 [2] CRAN (R 4.3.1)
 tidyr                  1.3.1      2024-01-24 [1] CRAN (R 4.2.3)
 tidyselect             1.2.1      2024-03-11 [2] CRAN (R 4.3.3)
 TMB                    1.9.11     2024-04-03 [1] CRAN (R 4.3.3)
 tradeSeq             * 1.14.0     2023-04-25 [1] Bioconductor
 TrajectoryUtils        1.6.0      2022-11-01 [1] Bioconductor
 tzdb                   0.4.0      2023-05-12 [2] CRAN (R 4.3.1)
 utf8                   1.2.4      2023-10-22 [2] CRAN (R 4.3.1)
 uwot                   0.2.2      2024-04-21 [2] CRAN (R 4.3.3)
 vctrs                  0.6.5      2023-12-01 [2] CRAN (R 4.3.1)
 viridis                0.6.5      2024-01-29 [2] CRAN (R 4.3.3)
 viridisLite            0.4.2      2023-05-02 [2] CRAN (R 4.3.1)
 vroom                  1.6.5      2023-12-05 [2] CRAN (R 4.3.3)
 withr                  3.0.0      2024-01-16 [1] CRAN (R 4.2.3)
 xfun                   0.44       2024-05-15 [1] CRAN (R 4.3.3)
 XML                    3.99-0.17  2024-06-25 [2] CRAN (R 4.3.3)
 xml2                   1.3.6      2023-12-04 [2] CRAN (R 4.3.3)
 xtable                 1.8-4      2019-04-21 [2] CRAN (R 4.3.1)
 XVector                0.40.0     2023-04-25 [2] Bioconductor
 yaml                   2.3.10     2024-07-26 [2] CRAN (R 4.3.3)
 zlibbioc               1.46.0     2023-04-25 [2] Bioconductor
 zoo                    1.8-12     2023-04-13 [2] CRAN (R 4.3.1)

 [1] /home/j.leary/r_packages_default
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

─ Python configuration ───────────────────────────────────────────────────────
 python:         /blue/rbacher/j.leary/py_envs/scLANE_env2/bin/python
 libpython:      /blue/rbacher/j.leary/py_envs/scLANE_env2/lib/libpython3.10.so
 pythonhome:     /blue/rbacher/j.leary/py_envs/scLANE_env2:/blue/rbacher/j.leary/py_envs/scLANE_env2
 version:        3.10.14 | packaged by conda-forge | (main, Mar 20 2024, 12:45:18) [GCC 12.3.0]
 numpy:          /home/j.leary/.local/lib/python3.10/site-packages/numpy
 numpy_version:  1.24.2
 
 NOTE: Python version was forced by use_python() function

──────────────────────────────────────────────────────────────────────────────